function [tau1,tau2,tau3] = directional_bary(V,T,x,y)
%
% calculate the barycentric of point (x,y) relative to all triangles
%
x1 = V(T(:,1),1); x2 = V(T(:,2),1); x3 = V(T(:,3),1);
y1 = V(T(:,1),2); y2 = V(T(:,2),2); y3 = V(T(:,3),2);

J  = (x2 - x1).*(y3 - y1) - (x3 - x1).*(y2 - y1);

% barycetric for (0,0);
J1 = (x2 - x ).*(y3 - y ) - (x3 - x ).*(y2 - y );
J2 = (x  - x1).*(y3 - y1) - (x3 - x1).*(y  - y1);
J3 = (x2 - x1).*(y  - y1) - (x  - x1).*(y2 - y1);

% barycetric for (0,0);
J10 = x2.*y3 - x3.*y2;
J20 = (x3 - x1).*y1 - (y3 - y1).*x1;
J30 = (y2 - y1).*x1 - (x2 - x1).*y1;

tau1 = (J1 - J10)./J;
tau2 = (J2 - J20)./J;
tau3 = (J3 - J30)./J;
end